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In this work we investigate the late-time stationary states of open quantum systems coupled to a 
thermal reservoir in the strong coupling regime. In general such systems do not necessarily relax to 
a Boltzmann distribution if the coupling to the thermal reservoir is non- vanishing or equivalently 
if the relaxation timescales are finite. Using a variety of non-equilibrium formalisms valid for non- 
Markovian processes, we show that starting from a product state of the closed system = system 
+ environment, with the environment in its thermal state, the open system which results from 
coarse graining the environment will evolve towards an equilibrium state at late-times. This state 
can be expressed as the reduced state of the closed system thermal state at the temperature of the 
environment. For a linear (harmonic) system and environment, which is exactly solvable, we are 
able to show in a rigorous way that all multi-time correlations of the open system evolve towards 
those of the closed system thermal state. Multi-time correlations are especially relevant in the non- 
Markovian regime, since they cannot be generated by the dynamics of the single-time correlations. 
For more general systems, which cannot be exactly solved, we are able to provide a general proof that 
all single-time correlations of the open system evolve to those of the closed system thermal state, 
to first order in the relaxation rates. For the special case of a zero-temperature reservoir, we are 
able to explicitly construct the reduced closed system thermal state in terms of the environmental 
correlations. 



I. INTRODUCTION 

Equilibrium states are typically discussed and derived 
in one of three settings or scenarios. In the more-common 
equilibrium (Gibbs) perspective, originally based upon 
classical ensemble theory, the entire system consisting 
of a system of interest plus its environment is taken to 
have some well-defined state or set of states, and upon 
coarse graining the environment, the system can appear 
thermal [1, 2]. In the less-common non-equilibrium per- 
spective, the environment is taken to be initially ther- 
mal, whereas the open system is allowed to dynamically 
relax from an arbitrary initial state into an equilibrium 
state [3-6]. This approach is referred to as the Langevin 
paradigm [7]. Both scenarios described above apply to 
situations where there is a clear distinction and separa- 
tion between the system and environment degrees of free- 
dom. When there is no clear distinction or the separation 
is not physically justifiable, as in a molecular gas where 
each particle is identical, a very different set of physical 
variables and different kind of coarse-graining measure 
need be considered. One can examine the behavior of the 
n-particle distribution functions and perform the coarse- 
graining (e.g., 'slaving' in [7]) on the Bogoliubov-Born- 
Green-Kirkwood-Yvon (BBGKY) hierarchy [8], This ap- 
proach is referred to as the Boltzmann paradigm. 

The equilibrium and non-equilibrium perspectives can 
be made to complement each other rather naturally 
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within the Langevin or open system paradigm. In the 
former case, Popescu et al. [ ] have shown that for an 
overwhelming majority of pure states of the system + 
environment (within a narrow energy interval), the re- 
duced density matrix is very close to the reduced density 
matrix corresponding to the microcanonical state of the 
system -I- environment (defined in the same energy inter- 
val). In their approach the comparison is done without 
explicitly determining an equilibrium state. The authors 
emphasize that for strong coupling, the equilibrium state 
is not of Boltzmann type, and yet their results are valid in 
this domain. It is important to note that dynamics does 
not play any role in their derivation; the entire argument 
is based on kinematics. The beauty of this approach is 
that one can explain the abundance of thermal-like states 
without referring to ensembles or time averages. 

Linden et al [ I] expands upon the approach of [1, 2| 
to demonstrate dynamical relaxation 1 under very weak 
assumptions. Specifically, they proved that any subsys- 
tem of a much larger quantum system will evolve to an 
approximately steady state. On the other hand Reimann 
[10] showed that the expectation value of any "realistic" 
quantum observable will relax to an approximately con- 
stant value. ([11] gave a clear analysis and unification of 
these two results.) Finally [12] proves relaxation over a 
finite amount of time both in the sense of [9] and [10]. 

Dynamical relaxation of an open quantum system has 
been studied in the limit of vanishing coupling to the 



1 See Sec. I A for the definition of the terms relaxation, equili- 
bration and thermalization as used in this work. There we also 
describe the meaning of the term equilibration as used in Refs. [9— 
12], which differ from our definition. 
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environment in [•!-<>]. In this limit the equilibrium state 
is shown to be of Boltzmann form in which case the result 
is called thermalization, rather than just relaxation. In 
our work reported here, we derive the equilibrium state of 
an open system coupled to a thermal reservoir explicitly, 
even in the strong coupling regime. Moreover for the N 
oscillator quantum Brownian motion (N-QBM) model we 
are able to show the relaxation of multi-time correlations 
of the open system as well. To do so we need to restrict 
the environment to be in a thermal initial state. 

Another difference between our work and [9-12] is in 
the methods and emphasis. We take the open quantum 
systems approach [6, 13-17] of assuming an environment 
(E) which the system (S) interacts with, keeping some 
coarse-grained information about the environment and 
accounting for its systematic influences on the system in 
a self-consistent manner. The time evolution of the open 
quantum system is in general governed by non-unitary 
dynamics. In contradistinction, [9, 10, 12] consider the 
unitary time evolution of the closed system (S + E) and 
then trace out the environment to get the system state. 
Both approaches are equally valid, each providing a dif- 
ferent perspective into the physics with different empha- 
sis. We will provide a more detailed comparison of our 
results to those in the literature in the discussion section 
at the end. 



A. Relaxation, Equilibration and Thermalization 

Before we present our approach, we want to define 
carefully what is meant by equilibration in this paper. 
To begin with let us consider a system in contact with 
two thermal reservoirs 2 at different temperatures. The 
system relaxes into a late-time steady state, which can be 
described by a reduced density matrix. All expectation 
values of system operators will also be time-independent 
at late-times. Yet there will be a steady heat flux from 
the hot reservoir to the cold reservoir through the system. 
This is an example of a non-equilibrium steady state. 

In general we define steady states via time independent 
density matrices: dp(t)/dt = and use the term relax- 
ation to describe the generic convergence of the reduced 
density matrix to a fixed but arbitrary state in the late- 
time limit. If the density matrix is diagonal in the energy 
eigenbasis of the system we call it a stationary state. An 
isolated stationary state is also a steady state, but this 
is not true for open systems with non- vanishing coupling 
to their environments. 

In this work we reserve the term equilibrium for sys- 
tems whose multi-time correlations can be derived from 
the thermal state of a possibly extended closed system 
which is governed by Hamiltonian dynamics. As a result 



2 We call an environment a reservoir if the environment has an 
infinite number of degrees of freedom, and a reservoir at constant 
temperature, a thermal reservoir. 



of our definition, equilibration implies relaxation but the 
reverse is not true. The thermal reservoir distinguishes it- 
self from other possible environments by the universality 
of its fluctuation-dissipation relation (FDR) 5 , detailed- 
balance conditions and Kubo-Martin-Schwinger (KMS) 
relations. In the vanishing coupling limit thermal reser- 
voirs lead to the thermalization of the system as defined 
below. However for non-vanishing coupling to a thermal 
reservoir the equilibrium state of the system does not 
need to be of the Boltzmann form 

*W = t£^*v (L1) 

The asymptotic states we derive in this paper in the 
strong coupling limit describe equilibration and not ther- 
malization. 

The term thermalization is reserved for the relaxation 
of the density matrix of a system to the Boltzmann form 
(1.1) irrespective of the initial state of the system. Ther- 
malization defined in this sense can take place only if the 
system-environment coupling is vanishingly weak. To be 
specific, one requires (1) decaying environmental correla- 
tion functions, defined in Sec. Ill, (2) an initially thermal 
reservoir and (3) vanishing relaxation rates 4 or, equiva- 
lently, vanishing environmental correlation functions. 




Environment 



FIG. 1. Depiction of a system embedded in its environment, 
with short-range interactions. The typical argument for ne- 
glecting the interaction energy is that in the macroscopic limit 
the boundary becomes immeasurable in relation to the bulk. 



3 As long as the environment is modelled after a physical system, 
fluctuations will be related to dissipation; hence there will be 
a FDR. However for general environments this relation depends 
on the specifics of the system-environment coupling. Thermal 
environments are unique in that the FDR does not depend on 
the details of the system and the coupling to the system [18]. 

4 To see a simple example of a relaxation rate consider the N-QBM 
model of Sec. II B for N=l. In the Markovian limit the damping 
kernel can be written as f(t, s) = 7oM<5(t — s), where 70 acts as 
the damping rate. 
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These conditions are customarily achieved by assum- 
ing short-range interactions and a relatively large system 
size, see Fig. (1). However this assumption is generally 
not justifiable for small systems as Fig. (2) suggests. In 
this paper, we address the stationary state of open quan- 
tum systems in contact with a thermal reservoir at tem- 
perature T = 1//3, without the assumption of a vanish- 
ing interaction strength and allow for finite relaxation 
timescales (). Relation (1.1) is known not to hold un- 
der these conditions [ ]. Phenomenologically, one can 
estimate the corrections we describe by the ratio of the 
relaxation rates 7 to the system's energy level splittings 
ft, or 7/ft. 5 

As thoroughly discussed in Ref. [19], this fact is of- 
ten overlooked in many circumstances, due to the effects 
of ancillary approximations such as the rotating-wave 
approximation, renormalization of environmentally- 
induced energy-level shifts and overly-simplistic models. 
As we explain in Appendix A, this fact may also be over- 
looked due to its absence in the case of classical, Gaussian 
noise. 




FIG. 2. Depiction of systems of decreasing particle number. 
For systems consisting of a small number particles, the argu- 
ment in Fig. 1 obviously does not apply. Furthermore, it is 
known that neglecting the interaction energy in these finite 
systems always results in infinite relaxation and thermaliza- 
tion times. 



described above. Here we refer to the result of these 
works using the terminology we defined above. 

B. Model and Assumptions 

We consider unitary dynamics of the closed system (C) 
described by the Hamiltonian He consisting of the sys- 
tem of interest (S) and its environment (E) with interac- 
tion (I) 

H c = H s + H E + Hi + H R , (1.2) 

where Hr contains all of the "renormalization" (R) ef- 
fects. The interaction generates environmental correla- 
tion functions, c.f. Eqs. (III. 4), (III. 8)), and we assume 
these correlations to be decaying functions. This assump- 
tion allows for irreversible dynamics in the open system. 
Implicit in this assumption is that the environment con- 
tains a continuum of modes (e.g. infinite volume). This 
latter assumption can be satisfied by coupling the system 
directly to field degrees of freedom that are uncountably 
infinite, such as the electromagnetic field. Note however 
that we do not assume the interaction Hamiltonian to be 
negligible compared to the system Hamiltonian. 

Finally, for mathematical simplicity we assume the ini- 
tial state of the system and environment to be uncorre- 
cted at t = 6 

e -/9H E 

Pc(0) = Ps(0)®T^y, (1-3) 

where the environment (a thermal reservoir) is in its iso- 
lated equilibrium state with partition function Z^{(3) = 
TrE[e _/3HE ], and the system (S) is in an arbitrary state. 

The assumption of a thermal state for the environment 
can be justified, for instance, by the approach of Popescu 
et al. [ ] in the weak-coupling limit, by giving the envi- 
ronment its own environment, without any restriction on 
the system-environment coupling strength. In this sense 
the work of Popescu et al., and those prior, serve as a 
pedagogical springboard for our analysis of strongly cou- 
pled systems. 



Finally, the term equilibrium is used in Ref. [9] to de- 
scribe what in our terminology are steady states and in 
Ref. [10] to describe what in our terminology are station- 
ary states. Both cases have been covered in Refs. [11, 12] 
with the single term equilibrium. These states do nec- 
essarily meet our more stringent criteria of equilibrium 



5 A well-known example is the density of states for an atom or 
molecule, which is necessarily interacting with the electromag- 
netic field to a degree which cannot be ignored when considering 
the Lamb shift, black-body radiation shifts, etc.. For optical 
frequencies, the emission rates of atoms are very small relative 
to their transition frequencies, and so these corrections are very 
small. However in other systems, such as condensates, these cor- 
rections can be of considerable size. 



C. Results 

It is well known that in the limit of vanishing inter- 
action strength, the open system relaxes to its thermal 
state [3, 6, 19, 22] 

e -/3H s 

lim i lim p s (t) = — — , (1.4) 

7->0t-Kx> ^s(P) 



6 The implication of initial correlations are considered in Ref. [20, 
21]: Correlated initial states are more physical, particularly in 
the early time evolution, but they have essentially no bearing 
on the mathematical results we derive herein, which are focused 
upon the asymptotic time evolution. 
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where p§(£) = Tte[pc(0] denotes the reduced density 
matrix and 7 a generic relaxation rate of the open system. 
Note that all relaxation rates are, at minimum, second 
order in the interaction, being primarily determined by 
the two-time correlations of the environment. 

In Ref. [ 23], it was shown to second-order in the inter- 
action, and for a single tensor-product coupling of system 
and environment operators, that the open system can be 
confirmed to relax to the reduced closed system thermal 
state 



lim p s (t) 

t— foe 



Tr F 



e -/3H c 



(1.5) 



We extend this proof to general system-environment cou- 
plings. For zero-temperature environments we demon- 
strate agreement with the ground state obtained from 
the time-independent Schrodinger equation. Moreover, 
we give a non-perturbative proof of Eq. (1.5) for the 
exactly-solvable model of iV-oscillator quantum Brown- 
ian motion (N-QBM), wherein the interacting system and 
environment are linear. In that model we are also able 
to rigorously prove that all multi-time correlations of the 
open system relax to those of the closed system thermal 
state with non-vanishing interaction. Correspondence of 
the multi-time correlations is an important consideration 
as, outside of the Markovian regime, the dynamics of the 
multi-time correlations cannot be generated by the dy- 
namics of the single-time correlations, as per the quan- 
tum regression theorem (QRT) [24]. 



The reduced, closed system thermal state 

It is important to emphasize that Eq. (1.5) pertains 
strictly to the open system S and not to the closed sys- 
tem (S + E), as equilibration requires not only a reservoir 
and late-time limit, but also a degree of coarse graining. 
As we show in Sec. II F, if one considers any individual 
mode of the environment, then its dependence upon the 
initial state of the system never decays. In this sense, 
information pertaining to the system's past is encoded in 
the environment, but only when considering the state of 
the closed system (S + E). Bowever, upon coarse grain- 
ing the environment by considering the time-evolution of 
a continuum of environment energies, and not one indi- 
vidual mode energy, then all dependence upon the initial 
state of the system is seen to decay away in time. In this 
sense, information pertaining to the system's past is only 
measurable for a finite span of time. 

The above statement is based on the fact that, while 
the open system experiences irreversible dynamics: dis- 
sipation, diffusion and decoherence, the closed system (S 
+ E) experiences reversible dynamics. Consider, for in- 
stance, the coupling of a mixed state of the system to a 
zero-temperature reservoir. Given unitary dynamics, the 
joint state of the system and environment cannot relax 
from a mixed state into a pure state (the ground state of 



the interacting theory). Bowever, the environment is ex- 
ceedingly large when compared to the system, and so the 
system's entropy, when spread out over every mode of the 
environment, can become immeasurable. This is a gen- 
eral phenomena of environmentally-induced irreversible 
dynamics: conserved quantities such as energy and en- 
tropy can flow into the environment, and owing to the 
overwhelmingly large number of degrees of freedom, be- 
come difficult to track or retrieve. 

The paper is organized as follows: In Sec. II we derive 
the equilibrium state for the linear N-QBM model. In 
Sec. Ill we extend our analysis to nonlinear systems via 
perturbation theory. In Sec. IV we summarize our results 
and compare them to relevant works in the literature 
and provide some new insights into the key issues. Some 
technical details have been provided and the notation is 
defined in the Appendices. 



II. LINEAR SYSTEMS 
A. The Lagrangian 

Our treatment of the N-QBM model is based on [25]. 
The model is that of a continuous and linear system 
with finite and countable degrees of freedom, with La- 
grangian L sys (X, X), bilinearly coupled, via a Lagrangian 
L; nt (X,x), to a linear environment with an infinite (and 
possibly continuous) number of degrees of freedom, with 
Lagrangian L onv (x,x). 

L = L sys (X, X) + L onv (x, x) + L int (X, x) + L rcn (X) , 



f / -. 



L = - (X T MX-X T CXJ + - (x T mx 
-x T g X + L rcn (X). 



(II.1) 



(II.2) 



We assume that the spring constant matrices C, c as well 
as the mass matrices M, m are real and positive definite, 
and can be considered in general to be symmetric. If nec- 
essary, one can relax the positivity condition and even 
consider time-dependent mass matrices, spring constant 
matrices and system environment coupling matrix g [2G] . 
Such a model environment can emulate any source of 
Gaussian noise with proper choice of coupling. To ensure 
that the free and interacting system are similar in behav- 
ior, we will also include the "renormalization" L rcn (X). 
Our choice of "renormalization" will be equivalent to in- 
serting the entire system-environment interaction in the 
square of the potential: 

L = - (X T MX-X T CX 

+ i (x T mx- [x-c _1 gX] T c [x-c-'gX]) , 

(II.3) 

since this keeps the phenomenological system-system 
couplings from changing. 
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B. The Langevin Equation 

For the linear system there are several formalisms 
which produce the same Langevin Equation. The most 
direct is via integrating out environment degrees of free- 
dom in the Heisenberg equations of motion [27] and then 
considering the symmetrized moments. Another is to 
consider the characteristic curves of the system + en- 
vironment's Fokker-Plank equation [26]. Finally, one 
can integrate out both the environment degrees of free- 
dom and the relative system coordinate A = X — X', 
while leaving only the average system coordinate X = 
(X + X')/2, in the double path integral of the reduced 
system propagator in the influence functional formalism 
[28] . In general (for nonlinear systems) there is no neces- 
sary correspondence between these formalisms and only 
the first may be well defined, but here the Langevin equa- 
tion is simply 

M X(t) + 2 [ ds-y(t, s) X(s) + C X(i) = £(f) - 2 -y{t) X . 
Jo 

(II.4) 

where 7 is the damping kernel and £ is the noise given 
by: 



7(M) 



T _ 1 cos(u;[i— si) _i 
2 -— y m 2 g 

2uS A 



£{t)=% (f(t)mx + f(i)p c 

, . _i sm(ut) _ 1 
{(t) = m 2 — - — -m 2 , 



u) = m 1 c m 



(II.5) 
(II.6) 

(II.7) 
(II.8) 



where f is the free Green's function of the reservoir posi- 
tions and us is the free reservoir frequencies upon diago- 
nalization. Note that the damping kernel is independent 
of the environment's initial state, whereas the properties 
of noise are determined by the environment's initial state. 

We consider the case in which the system and environ- 
ment are uncorrelated at t = and the environment is in 
its thermal state e~P HE /Ze((3). The noise has zero mean 
and the two time correlation is given by the noise kernel 



*(t,0 = (€(<) C(fT> 



(II.9) 



where the Gaussian average over the stochastic process 
£ is equivalent to tracing over the environment degrees 
of freedom. The noise and damping kernels satisfy then 
the fluctuation-dissipation relation (here in the Fourier 
domain) 



k(lu) = huj coth 



hid 



2k B T 

with the Fourier transform defined 



r+00 
/(w)= / dt, 



fit), 



(11.10) 
(11.11) 

(11.12) 



and where k is the (quantum) FDR kernel. Therefore, the 
problem is completely specified in terms of the damping 
kernel. 

Given that our damping kernel is stationary, the 
Langevin equation can be expressed in the Laplace do- 
main as 



[z 2 M + 2z 7 (z) + C] X(z) = [zMX 



■€(*), 

(11.13) 



where P = MX and (Xo,Po) correspond to the initial 
values at t = 0, and with the Laplace transform defined 



dte- zt 



/(*) 



(11.14) 



Formally, the solutions can be easily found by inversion: 

X(z) = G(z) [zM X + P ] + G(z) , (11.15) 
G(z) = [z 2 M + 2z7(z) + C] _1 . (11.16) 

Note that since our damping kernel is symmetric, i.e. 
■j(t, s) = ~f(t, s) T , the same will be true for the propaga- 
tor G(t, s) and its Laplace transform. It is also useful to 
consider the following representation: 



G(z) = M-5 



2z\(z) + n 2 



X(z) = Ms<y(z) M" 



: (11.17) 

(11.18) 
(11.19) 



where the eigenvalues of Q 2 coincide with the squared 
frequencies of the normal modes of the free system. Back 
in the time domain we have 

X(<) = G(t) M X + G(t) P + (G * £)(t) , (11.20) 

with * denoting the Laplace convolution, defined as 

(A*B)(t)= [ dsA(t-s)B(s). (11.21) 
Jo 

For more general Gaussian states, for which the sys- 
tem and environment are correlated, the noise can be 
correlated with (Xo,P<j) and the noise kernel modified. 
This is the case for the closed system thermal state given 
by the density matrix e _ ^ Hc /Zc(/3) which we investigate 
below. 



C. Single-time correlations in the closed system 
thermal state 



In this section we calculate the single-time correlations 
in the closed system thermal state of the N-QBM model. 
The partition function for the N-QBM model has been 
derived in App. C, Eq. (C.17). In the rest of the paper 
including the appendices we suppress the dependence of 
the paritition function on j3 for brevity of notation. As a 
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first step we take the logarithm of the partition function 
and write it as: 

logZ c = logZ E - ilrlogM- 1 - ^TrlogC 



differentiation of the partition function 
(xx T ) c 



TX 2ai0gZ C ,_ c -l g(xX T )cg T c -l 5 



dc 



^Trlog (m -1 ^^)" 1 ) + constant (11.22) ( X P T ) C = (p xT )c = °' 



r=l 



We begin by making a general observation. Consider 
the thermal state of a system described by a Hamil- 
tonian where the momenta appear only in the kinetic 
energy term of the form J^aPa/^ 171 - Then all correla- 
tions between position and momentum operators van- 
ish: (x a pi>) — 0. This can be seen by noting that 
all correlations are time-translation-invariant in equi- 
librium and forming the derivatives ^ (x a (t)xb(t)) and 
d ( t l f /) (xa(t)xb(t')) \ t _ t ,- This observation applies to N- 
QBM model. 

Let angular bracket with the subscript C denote ex- 
pectation values in the closed system thermal state. Ex- 
pectation values corresponding to the uncorrelated ini- 
tial state are denoted by attaching the subscript E to 
the bracket. For the purpose of partial differentiation, 
the partition function is to be regarded as a function 
of C, M, c, m, g and not (explicitly) of u). With a 
straight-forward application of Theorem 1, the reduced 
system correlations are given by: 



<XX T ) ( 



2d\ogZ c 



dC ' 

<xp t ) c = (px t ) = o, 
2aiogZc 



(PP T )c 



(3 dM- 



(11.23) 
(11.24) 
(11.25) 



The position-position and position-momentum correla- 
tions between system and reservoir modes are calculated 
similarly: 



(Xx T ) ( 



(XP T ) ( 
<P* T ) ( 



<*x t >; 

1 SlogZc 

(PX T )" = 
(xP T )c 



<XX T ) ( 



T -1 
g C 



= 0, 
0. 



(11.26) 

(11.27) 
(11.28) 



To calculate the momentum-momentum correlations be- 
tween system and environment we take the time deriva- 
tive of (X(t)p T (i)) c and set it to zero. Since in the 
closed system thermal state all expectation values are 
time-independent, we know that there is in fact no de- 
pendence on time. Using the equations of motion it is 
straight-forward to show that: 

(Pp T ) c =M (Xx T ) c c-M (XX T ) C g T . (11.29) 

The environment correlations can be calculated by direct 



, T > 2 dlogZc 

\PP ) C £ 



(11.30) 
(11.31) 

(11.32) 



Now we are in a position to determine all the single- 
time correlations of the interacting theory in the closed 
system thermal state. Since the equilibrium state 
is stationary these single-time correlations are time- 
independent. The details for some of these formulae are 
provided in App D. All the nonzero correlations are given 
by: 



(xx T ) c = iGM + |£^) 



1 



( 

2 



<PP T ) C = - M-z/ 2 MG(^ )M 



(11.33) 
(11.34) 



P 



jr(M-^MG(iy r )M 



P r=l 

(11.35) 

(Pp T ) c =M(Xx T ) c c-M(XX T ) c g T , (11.36) 

(PoPo) c = <PoPo)e ( IL3? ) 

2 °° 

P r=l 

(x xT) c = (x x£) E + c^g (X X) c g^ 1 (11.38) 
2 



P r=l 



- B E ^ (f (*) g 6("r) g T f (^r) m C- 1 ) 
P r=l 

2 °° / 

- ( c " lmf K)gGWg T %)mc 

P r=l 

where f is the Laplace transform of the free reservoir 
propagator given by Eq. (II. 7) and v T = 2-irr/HP are the 
Matsubara frequencies. 



D. Equivalence of single-time correlations for the 
open system 

In this subsection we show that the single-time correla- 
tions of system variables for the uncorrelated initial state 
are asymptotically identical to the single-time correla- 
tions corresponding to the closed system thermal state. 
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We start by calculating the variances for the closed sys- 
tem thermal state. The requirement that G(t) is a decay- 
ing function means that the Laplace transform G(z) is 
analytic in the right half-plane. Hence G(— iuj) is analytic 
in the upper-half plane. On the other hand coth(/3fiw/2) 
has simple poles on the imaginary axis at the Matzubara 
frequencies v^. The summations over r in Eq. (11.33) can 
be written as a contour integral using Cauchy's theorem: 



<XX T ) ( 



pn/2 

— — x ■ 
2m /? J c 



dzcoth(/3hz/2)G(-iz) 



(11.39) 



A similar argument can be used to derive: 



(PP T )c 



h 

2^ 



duui coth(/3fiw/2)Im G(-itu 



(11.41) 



Eqs. (11.40,11.41) are identical to the results obtained by 
[ ] for the asymptotic values of variances corresponding 
to an uncorrelated initial state. Therefore we have proven 
that the single-time correlations of the open system relax 
to those of the closed system thermal state. 



E. Equivalence of multi-time correlations 



The contour of integration is chosen to encircle the upper- 
half plane in a counter-clockwise direction. The poles on 
the imaginary axis at Matzubara frequencies v r for r > 1 
are encircled, but only half of the pole at the origin is 
enclosed. The arc of the contour does not contribute to 
the integral when the radius is taken to infinity. Hence 
we can write this expression as an integral on the real 
line. Furthermore, by the symmetry of the integrand, 
the real part vanishes and the integral is given by: 



2- 



+ 00 



dwcoth(/3fio;/2)Im G(-VjS) . 

(11.40) 



In this section we generalize the results of the pre- 
vious section to include multi-time correlations. We 
begin by calculating the two-time correlation function 
(X(t) X(i') T ) c using the trajectories obtained from the 
Langevin equation. Note that for the closed system ther- 
mal state this quantity is stationary. To simplify the 
proof we make use of this observation and take the late- 
time limit of the closed system thermal state as well with- 
out loss of generality. This trick makes the comparison 
of the two cases easier and reduces the amount of com- 
putation. 

The dynamics of the system is given by the solution 
(11.20) of the Langevin equation which is valid for any 
initial state. The dependence on initial state is hidden in 
the correlations between Xo, Po and £(t). The two-time 
position correlation is given by 



<X(t) X(0 T ) C = G(t) M (X X r ) c M G(i') + G(t) <P P?) c G T (t') 



G(t) M / ds' (X £( s ') T ) c G(t'-s')+ fd s G(f- s )((( s )Xj) c MG(t') 
Jo Jo 

G(i) fds> (P ^s') T ) c G(<'-,s')+ fdsG{t- S ){Z{ S )Pl) c G{t') 
Jo Jo 



f*da fds' G(t-s) t(s') T ) c G(t'-s') 
Jo Jo 



(11.42) 



As mentioned earlier unlike the uncorrelated initial state 
the terms in the second and third lines do not vanish in 
the closed system thermal state. We consider the case 
where 



lim G(t) 

t—>oo 



lim 7(i) = . 

t— ¥00 



(11.43) 



This is the criteria for dissipative dynamics. Under these 
assumptions the first two terms in Eq. (11.42) vanish in 
the late-time limit for any initial state. The terms in the 
second and third lines have one factor of G(t) or G(t) 
that goes to zero in the late-time limit multiplied by a 



convolution integral. In App. D we show that these con- 
volution integrals are finite. Hence the terms in second 
and third lines also vanish asymptotically. Finally we 
show the equivalence of the term in the last line for the 
uncorrelated and thermal initial states at late times in 
App. E. 

The comparison of more general multi-time correla- 
tions can be done similarly using the trajectories of the 
Langevin equation. The above example demonstrates 
how in the late-time limit the effects of initial conditions 
of the system die out and the noise statistics of both 
preparations converge. The equivalence at the level of 



trajectories ensures that all the multi-time correlations 
will be identical. 

Let us reiterate the result we just obtained: a lin- 
ear system linearly coupled to a linear thermal reservoir 
(with uncountably many degrees of freedom) at inverse 
temperature (3 does relax to the equilibrium state de- 
scribed by (1.5). This state is different from the Boltz- 
mann state given by (1.4) whenever the interaction be- 
tween the system and environment is not negligible. 
Moreover the multi-time correlations of system observ- 
ables also relax to their corresponding values in the closed 
system thermal state. 



F. The effect of coarse graining 

Up until this point we only focused on the system de- 
grees of freedom. Now we turn our attention to the envi- 
ronment. Following Ref. [25, 26], the trajectories of the 
environment oscillators, as driven by the system oscilla- 
tors, are given by 

x(i)= [f(()mx(0) + f(t)p(0)l +f(t)*gX(t), (11.44) 

in terms of their free propagator f(i) and frequency ma- 
trix u> given by Eqs. (11.7,11.8). Into Eq. (11.44) we sub- 
stitute the system trajectories, which are damped oscil- 
lations driven by noise for the continuum environment: 



the system is given by 

/ x f d P(0) 1 f d 2 „ d 



h k (t) 



9kfk(t) 



(u,£-ffT + 47o 2 



2, ,2 



h k {t) 

(11.49) 
(11.50) 



plus terms that decay exponentially and the terms which 
depend upon the initial state of the environment. The 
function hk(t) oscillates forever, the same as /fe(t), and 
therefore the environment retains information pertaining 
to the initial state of the system forever. However, this 
information is not measurable forever. The system only 
interacts with the integrated trajectories, which resolve 
to a convolution of the damping kernel and open-system 
propagators. 



g T x(t) = -2 j{t) * G(i) M X(0) + G(t) P(0) 



(11.51) 



and upon integrating over a continuum of environment 
frequencies (here performed by multiplication with the 
infinite matrix g T ) the oscillatory terms decay in time. 
Thus the late-time limit and coarse graining together are 
responsible for the erasure of all information pertaining 
to the initial state of the system. 



X(t)= G(i)MX(0) + G(i)P(0) +G(i) *€(<)• 

(11.45) 

We then find the environmental dependence upon the 
initial state of the system to be 



x(f)=f(t)g* G(i)MX(0)+G(t)P(0) 



(11.46) 



III. GENERAL SYSTEMS 

A. Steady state 

The time-evolution of the reduced density matrix of the 
open system can be generated by a perturbative master 
equation 



Ps(t) = £(*Hps(*)}> 



(in i) 



with all additional terms only dependent upon the initial 
state of the environment. The system-dependent terms 
correspond to a convolution of harmonic oscillations of 
the environment with non-locally damped oscillations of 
the system. Resolving these integrals leads to some terms 
which oscillate with environment frequencies uj and do 
not decay. 

As a simple example, consider the local damping of a 
single system oscillator. The open-system propagator or 
Green's function is given by 



n 



Am 
ft 2 



(11.47) 
(11.48) 



The environment's dependence upon the initial state of 



where the Liouville operator can be expanded in terms 
of the interaction Hamiltonian by a variety of methods 

[22, 29-31]. 



C(t) = Co + d(t) 
A>{p} = -i[H s ,p] , 



c 2 (t) + 



(III.2) 
(III.3) 



In general, C\(t) can be absorbed into the system Hamil- 
tonian and so we will primarily concern ourselves with the 
second-order term. For simplicity we will assume there is 
no degeneracy or near-degeneracy in the system energy 
spectrum. 

Expanding the interaction Hamiltonian in terms of sys- 
tem L n and environment 1„ operators 



H, 



L 



(III.4) 
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the multivariate master equation can be represented [22] 
£2 P = ^2 [L«, P (A nm o L m ) t - (A nm L m ) p] , 

nm 

( IIL5 ) 

where the A operators and product define the second- 
order operators 

(A nm oL m )(i)= / dsa nm {t,s) {Go(t, s) L m (s)} , 
Jo 

(III.6) 

in terms of the zeroth-order (state) propagator of the 
system 



9o(t,s){p} = e 



- „-t(t-s)H E 



pe 



+i(t-s)H s 



(III.7) 



and the (multivariate) environmental correlation func- 
tion 



l (t,s) = (l n (t)l m (s)) I 



(III.; 



The second-order operator can be expressed as the 
Hadamard product 

A nm o L m \uii>) = A(uj t ~uJi') (u)i\ L TO |w,v) , (III. 9) 

and, in the late-time limit, the second-order coefficients 
resolve 



1 



4m(w) = -5w(w) - lV 



*a nm (u), (111.10) 



where dc(u>) denotes the Fourier transform of the station- 
ary environment correlation function a(t— s) = a(t, s), V 
the Cauchy principal value and * the appropriate Fourier 
convolution. 

With the multivariate master equation detailed, we can 
prove relation (1.5) to second order in the interaction. 
This generalizes the univariate proof in Ref. [23], which 
considered a single tensor-product interaction between 
the system and environment. As the proof is straightfor- 
ward in either case, we will give an outline and focus upon 
differences which arise in the multivariate treatment. 

We are looking for the stationary state pp, such that 



C{pA = 0, 



(111.11) 



we know from detailed balance that the zeroth-order 
stationary state is the thermal state (1.4), e.g. see 
[22]. Second-order corrections can be generated from the 
second-order master equation via canonical perturbation 
theory. More explicitly, we have 



(^|£ 2 { e ^ H }|^- 



U>i — Ldj 

(111.12) 

but only for the denoted off-diagonal perturbative cor- 
rections (in the energy basis |w)). As explained in 
Ref. [32] , due to unavoidable degeneracy, specifically that 
the diagonal elements are all stationary to zeroth-order, 
the second-order master equation cannot determine the 



second-order corrections to the diagonal elements of the 
density matrix. Fortunately this does not greatly hamper 
our proof of correspondence. By a simple application of 
the multivariate master equation to Eq. (III. 12), we eas- 
ily obtain these second-order corrections to the thermal 
state of the system. 



B. Equilibrium state 

We wish to compare the straightforward expansion of 
(III. 12) to the reduced closed system thermal state at 
second order, and so we require a perturbative expansion 
of (1.5). There exists such a perturbative expansion of 
exponential matrices utilizing the identity 



de 



-e A+£B = e 



A+eB / d U e-"( A + £B »B e + , '( A+£B ), 
Jo 



(111.13) 

to obtain an operator- Taylor series in the perturbation 
eB. After a fair amount of simplification, one can deter- 
mine the second-order stationary state to be 



-/3H S 



+ e" 



-PHs 



(III. 14) 

' /V /"V'(H I (-0 / )H I (-z^")) E , 
Jo Jo 

in terms of the complex-time operators 

HtC-ijS) = e+^Me) H x e-^H+H*) , (111.15) 

where the noise average is taken with respect to the free 
thermal state of the environment and factors inside the 
environmental trace have been written to suggest their 
correspondence with the environmental correlation func- 
tion evaluated at imaginary times. 

The double integrals in Eq. (III. 14) reduce to 

f V f dP" a nm (-iP', -iP") L n (-ip') L m (-iP") , 

(111.16) 

in terms of the complex-time operators 

Lt-^Ee+^Te-" 115 . (111.17) 

After a Fourier expansion of the complex-time correla- 
tion functions, expressions (III. 12) and (III. 16) can be 
compared term-by-term in the energy basis wherein the 
imaginary-time integrals of Eq. (III. 16) can be resolved as 
the master equation operators were. Though the two ex- 
pressions will then be composed of the same objects, they 
will not immediately appear to be equivalent. The final 
step is to apply the relevant multivariate Kubo-Martin- 
Schwinger (KMS) relations (also found in [22]) 



a(+cu) = a T (— uj) e ?=ct*(— u)e ? 



(111.18) 



and then one can see that the two expressions are equiv- 
alent, even in their missing diagonal perturbations. 
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As previously discussed, Eq. (III. 12) is missing diag- 
onal perturbations due to degeneracies inherent in all 
perturbative master equations. The same discrepancy 
in Eq. (III. 16) stems from the perturbative expansion of 
the equilibrium state being inherently secular in (3. Both 
expressions are equivalent and quite conveniently they 
both lack the same second-order corrections to the di- 
agonal entries of the density matrix. Therefore, as far 
as the second-order dynamics is concerned, our proof of 
correspondence is complete. 



C. Zero- Temperature Analysis 

Though correspondence was established, the previous 
analysis was seen to be insufficient for calculation of low- 
temperature equilibrium states of the open system. How- 
ever, as we shall now show, at least for zero-temperature 
noise, it is still possible to easily construct the reduced 
closed system thermal states in terms of the same en- 
vironmental correlation functions which occurred in the 
previous analysis. The following relations were applied 
towards the inspection of two-level atoms interacting via 
a zero-temperature quantum field in [33]. 

In the zero-temperature regime we can apply mundane 
perturbation theory to derive the stationary-state per- 
turbations. One merely considers the perturbed ground 
state of the system + environment 



= *l> + + V> 2 + 
0o = |0) ® |0) E , 

and then traces out the environment 



(111.19) 
(111.20) 



Pp = 



+ (0 2 0S + 1 t 1 +0 o t 2 ) E + --- , (111.21) 



where we neglect the first moment of the reservoir as 
previously discussed. Without loss of generality let us 
set the ground-state energy of the system to zero. The 
second-order reduced ground state of ordinary perturba- 
tion theory provides the following stationary state per- 
turbation 



nmk 



L m bk) (oj k \ L„ \u 



j i > 



(111.22) 

where Zo(J3) is the partition function of the free system 
and with the off-diagonal (and non-resonant) coefficients 
given by 



-An 



~^A mn (uj kl ) - e-^iA mn {w kj ) 



(in.23) 



+ An 
where oj, 

master equation coefficients in (III. 10). 'An' denotes the 



ujj and A nm (uj) are the second-order 



anti-Hermitian part; the Hermitian and anti-Hermitian 
parts are defined 



He[Q 
An[Q 



inn} — 2 {Qnm "i" Q m 



mn} — 2 (Qnm Qmn) 



(111.24) 
(111.25) 



and for univariate noise (one collective coupling to the 
reservoir) the Hermitian and anti-Hermitian parts are 
simply the real and imaginary parts. In either case the 
anti-Hermitian part of (III. 10) is the second term. These 
off-diagonal perturbations perfectly correspond to the re- 
sults of the previous section upon introducing the appro- 
priate Boltzmann weights as we have. Moreover we also 
have the diagonal (and resonant) coefficients determined 
to be 



ijk \uii—Ui 



(111.26) 



An 



e fiuk A nm (u!ik) + e A rnn (u>ki) 



duji 



duii 



where here the Boltzmann weights are guessed, as 
these relations have only been derived here for zero- 
temperature. Therefore Eq. (III. 26) is exact for 
zero-temperature and our best guess for the positive- 
temperature coefficients: it has the correct functional de- 
pendence upon the Boltzmann weight and fourth-order 
master equation coefficients. At worst this is an interpo- 
lation of the zero and high-temperature states. 



IV. DISCUSSION 

In this work we investigate the equilibrium states of 
open quantum systems from dynamics / non-equilibrium 
point of view. We show that starting from a product state 
(1.3) the open system which results from coarse graining 
the environment will evolve to a late-time steady state. 
This state can be expressed as the reduced state of the 
closed system thermal state at the temperature of the en- 
vironment, i.e. Eq. (1.5). This result is important when 
the system-environment coupling is not negligible 7 , or al- 
ternatively, when relaxation rates are not insignificant in 
relation to the system frequencies. In this case the sta- 
tionary state of the system (1.5) differs from the canonical 
Boltzmann state (1.4). 8 One might argue that this state 
is the closest one can get to thermalization in the strong 
coupling regime. 9 However in this paper we use the term 



7 Based on the discussion of Fig. 2, we expect our results to be 
most relevant to small systems. 

8 In this paper we have not focused on the nature of this difference. 
A quantification in terms of the Hamiltonian of mean force for 
the special case of an Ohmic environment is given by Hilt et al. 
[34]. We intend to address this issue in our future work. 

9 Alternatively one could define this state to be the thermal state 
in the strong coupling regime. However this state depends on the 
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equilibrium state for Eq. (1.5) and reserve the term ther- 
mal state to the standard Boltzmann form (1.4). 

Our proof is exact for the linear model and to sec- 
ond order in interaction strength for nonlinear models. 
Moreover, for the exactly solvable linear case we prove 
the equivalence of multi-time correlations. The issue 
of multi-time correlations in the context of equilibra- 
tion/thermalization seems to be mostly ignored in the 
literature. We argue that multi-time correlations are im- 
portant outside the Markovian regime, as was pointed 
out in [28]. For instance, the relaxation of multi-time 
correlations cannot be deduced from the relaxation of the 
reduced density matrix of the system, neither can the ex- 
plicit value of the multi-time correlations be derived from 
the equilibrium state, if the dynamics is non-Markovian. 
In this respect our analysis of the linear N-QBM model 
provides insight into equilibration phenomena beyond the 
density matrix formalism. 

A complete proof, which would be non-perturbative 
for non-linear systems, would have to be very different 
than the second-order proof presented here. Our nonlin- 
ear proof, though very general in its application to dif- 
ferent systems and environments, is not robust enough 
for non-perturbative multi-time correlations. It is not 
immediately clear how such a proof could be attempted, 
whereas the elegance of the final result makes the possi- 
bility of its existence seem reasonable. 

An analogous proof for classical systems should be 
attempted by coarse graining the symplectomorphic 
(Hamiltonian) time evolution of the system and envi- 
ronment in much the same way that quantum master 
equations result from coarse graining the unitary time 
evolution of the system and environment. Unfortunately 
the literature on such an analog is not well developed 
(e.g., it would involve higher-order Fokker-Plank equa- 
tions which might only perturbatively preserve probabil- 
ity) and this would be more mathematically challenging 
than the quantum proof. Note that the h — > limit of the 
quantum results obtained in this paper yield the corre- 
sponding classical results, as has been argued in App. A. 

An essential ingredient of our proofs is the continuum 
limit for the environment. For a finite environment the 
t — > co limit of the reduced state does not exist within 
the formalism presented here and another ingredient is 
necessary to ensure relaxation to equilibrium. Having 
classical molecular dynamics in mind, we entertain the 
possibility that quantum chaos might be one avenue to 
explore. 

On the other hand we can consider a large but finite 
environment. It can be argued that for any relevant times 
t > the effect of an infinite reservoir can be approxi- 
mated arbitrarily closely by a large but finite reservoir. 
Then equilibration is observed for the time-interval be- 



specifics of the reservoir and the coupling to the reservoir. Hence 
it is not specified by the system parameters alone and referring 
to it as the thermal state is, in our opinion, misleading. 



tween the relaxation time and the recurrence time. Note 
that this interval is huge for a large environment, since 
the recurrence time grows very rapidly with the number 
of degrees of freedom. As a result the system stays close 
to its equilibrium most of the time. This interpretation 
helps us touch base with the results of [9, 10, 12] where 
relaxation in finite systems is proven for time averaged 
quantities. 

A. Comparison with recent literature 

To put this work in developmental context, here we 
compare more specifically our results to that of Linden 
et al. [9], Reimann [10], and Short and Ferrelly [12] 10 . All 
these works have in common with us the set-up of a small 
system coupled to a large environment and relaxation is 
achieved dynamically via time-evolution. A major dif- 
ference is the choice of initial conditions: they allow for 
any initial state, which is spread over sufficiently many 
energies, whereas we restrict our environment to be in 
a thermal state. In turn we can derive the form of the 
equilibrium state explicitly. 

Unlike what is done here these authors all make the as- 
sumption of non-degenerate energy gaps (this assumption 
is relaxed to a certain degree in [12]) and assume finite 
dimensional Hilbert spaces. The linear model we solved 
exactly here has infinitely degenerate energy gaps and we 
considered a reservoir consisting of an infinite number of 
degrees of freedom. Ref. [9] considers only pure states 
for the closed system (in the sprit of [1, 2]). Finally they 
all define relaxation in terms of time averaged quantities, 
i.e. systems behave as if they are in their steady state 
most of the time. Ref. [ ] also provides upper limit for 
the relaxation time. 

The proofs of [9-12] rely on the much greater dimen- 
sionality of the Hilbert space of the environment com- 
pared to that of the system. The system + environment 
state is propagated as a whole using unitary dynamics. 
The fact that the environment is large is utilized in the 
tracing out of the environment at the end of time evolu- 
tion. In this derivation the effect of the environment on 
the system dynamics is not so easily accessible. 

In our proof, the fact that the environment consists of 
a large number of degrees of freedom manifests itself in 
the form of its decaying correlations. These correlations 
in turn determine the non-unitary aspects of the open 
system dynamics. We use this non-unitary open system 
dynamics to evolve the reduced state of the system to its 
equilibrium state. In particular we do not refer to the 
state of the closed system explicitly 11 . Our derivation is 



See Sec. I A for the clarification of the different use of the term 
equilibration in the literature and here. 

Except for Sec. II F, where we do look at the individual environ- 
mental modes just to make the point that the closed system (S 
+ E) does not equilibrate. 
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more in the idioms of open quantum systems paradigm, 
where the influence of the environment on the system 
dynamics can be continuously monitored and explicitly 
expressed (e.g., consistent backreaction from the environ- 
ment is fully embodied in the influence functional [ ]). 

Relaxation is demonstrated in [9, 10, 12] for very gen- 
eral Hamiltonians, including strong coupling between the 
system and the environment. In their derivation the 
strong coupling regime does not present any extra dif- 
ficulty. In the open system approach we adopted in this 
paper strong coupling is difficult to handle. On the other 
hand, as a benefit of our method we can describe the 
nature of the equilibrium state, i.e. Eq. (1.5), besides 
proving its existence and uniqueness. 



Appendix A: The triviality of classical, Gaussian 



While Ref. [19] gives many cases in quantum mechan- 
ics in which the effect of system-environment coupling on 
the equilibrium state may be overlooked, here we would 
like to motivate the fact that this point is often over- 
looked in the classical regime as well, perhaps due to the 
ubiquitous employment of Gaussian noise. Let us con- 
sider the Hamiltonian of a system coupled linearly, via 
the system operator L, to an environment of harmonic 
oscillators, indexed by k, which mock our Gaussian noise 
[13, 35]. 



H c = H s + ]T 



H s 



E 

k 



2m k 

jpL 

2m k 



m k 0J% 2 



m k u 2 k 



L E gkXk 

k 

m k u\ 



H, 



(A.l) 



where the linear interaction is included in the square of 
the environment potential as a means of "renormaliza- 
tion" . Otherwise, the influence of the environment effec- 
tively introduces a negative L 2 term proportional to the 
cutoff into the system Hamiltonian when considering the 
open-system dynamics. 

Tracing over the environmental degrees of freedom is 
equivalent to integrating over the environmental dimen- 
sions in phase space, 



Tr, 



n fa / 



dpk ■ 



(A.2) 



where classically-speaking, x k and p k are independent, 
commuting variables. Therefore, in the classical and 
Gaussian model, relations (1.4) and (1.5) are equivalent 
as tracing over the environmental degrees of freedom con- 
stitutes a trivial Gaussian integral in phase space. The 
classical result can also be reached as the H — >• limit of 
the quantum result. This limit is most straight-forward 



when applied to the Wigner function [ ] defined as: 
1 



W(x,p) = ldue* pu p{x-u/2,x+u/2) . (A.3) 
2i:h J 

The description in terms of the Wigner function is equiv- 
alent to the density matrix approach. Hence the Wigner 
function contains complete information about the quan- 
tum system. As a result the Wigner function should not 
be treated as a phase space distribution, since it can as- 
sume negative values. However the H — ► limit of the 
thermal state Wigner function is well-defined and gives 
the classical Boltzmann distribution function: 



lim Wb(x,p) 



(A.4) 



For classical open systems it is well known that if 
the system + environment is in a thermal state of the 
full Hamiltonian, which includes the system-environment 
coupling, then the reduced distribution of the system is in 
general not the thermal distribution of the system Hamil- 
tonian alone. The term potential of mean force is used in 
chemical-physics literature for the quantity that replaces 
the Hamiltonian in the familiar Boltzmann distribution 
[ ]. The linear reservoir is a special case where the po- 
tential of mean force coincides with the system Hamilto- 
nian. The potential of mean force is defined by 12 : 



H*(X,P) 



1 U k IdxkJdp k e-^^ p ^ 
> ° g EL' /da* /dp*/ e-/> H *<*.p) 



(A.5) 



To the best of our knowledge, the asymptotic time evo- 
lution of a general classical open system, with a nonlinear 
environment initially in its thermal state, is not known. 
We conjecture that the reduced system is asymptotically 
described by e _ ^ H as described in the previous para- 
graph, and as would follow from (1.5). In this paper we 
provide a proof of the analogous statement for quantum 
systems to second order in interaction strength. Obvi- 
ously, our second-order proof extends to classical systems 
which can arise in the limit h — > 0. For linear systems we 
have an exact proof, and unlike its classical counterpart, 
the quantum linear case is highly nontrivial. 



Appendix B: Theorems on matrix derivatives 

Notation and Remarks: A letter in bold like A indi- 
cates a matrix. Referring to an element of the matrix we 
use subscripts: A a b. The inverse of the matrix is indi- 
cated by A -1 . An element of the inverse matrix is writ- 
ten as (A _1 ) a b to avoid confusion with 1/A a b. Transpose 



12 In most treatments Hg is absent. In that case even for linear 
reservoir H* differs from Hg by a frequency "renormalization" . 
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of the matrix is denoted by A T . Tr without a subscript 
indicates ordinary matrix trace. Trc indicates quantum 
mechanical trace over the closed system Hilbert space. A 
systematic study of matrix derivatives including some of 
the theorems below is given by [ ] . 

Before proceeding to the derivations we clarify a math- 
ematical subtlety. The theorems derived in this appendix 
will mostly be applied to symmetric matrices for which 
A a b = Aba- When taking the derivative of such a matrix 
with respect to one of its elements one can adopt two dif- 
ferent conventions. If the derivative is taken under the 
constraint that only symmetric variations of the matrix 
is allowed the result is: 



9A 
dA, 



ab 
cd 



<5ac^bd + <$ad<5bc(l - ^ab) 



(B.l) 



On the other hand if independent variations of all ma- 
trix elements are allowed the second term in the above 
equation is absent. In the following theorems we adopt 
the second convention. 

Theorem 1. Consider a system in a thermal state at 
inverse temperature f3 described by a Hamiltonian with 
parametric dependence on a set of variables {A„}. Then 
the expectation value of the derivative of the Hamiltonian 
with respect to these parameters can be calculated from 
the partition function by: 



OH 

8X n 



Tr ( 



m e-? H 



l_d_ 



ln(Z) . (B.2) 



Proof. In this proof we will make use of the following 
operator identity valid for an arbitrary operator O: 



d 
dX, 



dO 



e° = / due 110 — 



(B.3) 



/o d\ n 
Using this formula we can write the RHS of Eq. (B.2) as: 



pz 



:Trr 



8X n 
■l 



du e 



(B.4) 

-U0H d/j-H (i_„)^h 







BXr, 



We use the cyclic property of trace to get: 

r-l 



1 



:Tr c 

OH 

dX 



o 



oX n 



riH 



,-/8H 



(B.5) 

(B.6) 
(B.7) 
(B.8) 



Theorem 2. For a matrix A 

TrlogA = logdet A. (B.9) 

Proof. Trace operation is basis-independent. In the basis 
in which A is diagonal log A is also a diagonal matrix 
with entries log a n where a n are the eigenvalues of A. 
Taking the trace gives: 



Tr log A = lo S On = log a n J 

n \ n / 



(B.10) 



The last expression is recognized to be log det A since the 
product of eigenvalues equals the determinant. □ 

Theorem 3. For an arbitrary number of matrices A^ 
indexed by k, the following is true: 



Trlo g mA fe ) =^Trlog(A fe ) 

V k J k 



(B.ll) 



Proof. To show this equality we make use of Theorem 2, 
the well known fact that the determinant of the product 
of matrices equals the product of the determinants and 
properties of ordinary logarithms: 



Tr log (j[ A fc ^J = log det (j[ A k ^j , 
= log ^JJ det A^j , 



^2 log det A k , 

k 

5>rlog(A fe ). 



(B.12) 

(B.13) 
(B.14) 
(B.15) 

□ 



A corollary of this theorem is the fact that Trlog is 
invariant under any permutation of its arguments. 

Theorem 4. Consider a matrix A and a parameter X. 
Then: 



0_ 

OX 



Tr log A = Tr 



_ x dA 
~8X 



In particular: 



8A 



TrlogA = (A _i ) 



1\T 



(B.16) 



(B.17) 



where -J^ is defined as the matrix obtained by differenti- 
ating with respect to the entries of matrix A. 



n i c 



Proof. Let A = 1 + B and use 

log(l + B) = B - B 2 /2 + B 3 /3 ■ 



(B.18) 



□ to write the LHS of Eq. (B.16) as: 



14 



dX 



Tr 



B 



B 2 

2 



B 3 

3 



Tr 



9A "2 ( 9A B + B 9A ) + 3 ( 9A B 



B aA B + B V 



(B.19) 



Using the cyclic property of trace we obtain: 



Tr 



(l-B + B 2 -B 3 



(B.20) 



Note that dB/dX = dA/dX and 



1 — B + B — B 



= (1 + B)- 1 = A" 1 , (B.21) 



which proves Eq. (B.16). To prove Eq. (B.17) let A 
A ab . 



8 



dA 



- Trios; A = Tr 



ab 



_! dA 



<9A ab _ 
dA dc 



-E( A ^)cd ( ; Aa!i - 



cd 



(B.22) 
(B.23) 
(B.24) 



□ 



Theorem 5. Let A be an invertible matrix and X a pa- 
rameter. Then: 



dJ^_ dA 
OX d\ 



In particular: 

d(A-% h 
OA 



-(A- 1 ) am (A- 1 ) nb . 



(B.25) 



(B.26) 



Proof. We write A^ 1 = A -1 A A -1 , and differentiate 
both sides with respect to A. Looking at an element of 
this matrix equation we have: 



ox 



E 



dX 



cd 

dX 



'■A cd {A- x ) dh + (A- 1 ) 



dA cc 
' dX 

d(A-% h 
dX 



(A- 1 ) db + (A- 1 ) ac A c 



d(A~ 



OX 



(B.27) 
(B.28) 



This proves Eq. (B.25). For the proof of Eq. (B.26) we 
set A = A mn . □ 

A corollary of this theorem is the following identity 
valid for independent matrices A^: 



d 
&Al 



Tr [A 2 Af 1 A 3 ] 



(A^AaAaAT/ 1 ) 



(B.29) 



Appendix C: N-QBM Partition Function 

In this section we calculate the partition function of 
the N-QBM model. Our treatment mimics and general- 
izes that of Weiss [16], which treats one system oscillator 
only and does not allow for interactions among reservoir 
oscillators and non-diagonal mass matrix. 13 The parti- 



tion function has an imaginary-time path integral repre- 
sentation given by: 



Z c = ^DxDXcxp (-S (E) [x,X]/ft 



dri( E )(r), 



i (E) (r) = i (X T MX + X T CX 



(CI) 
(C.2) 
(C3) 

+ i(x T mx+[x-c- 1 g X] T c[x-c- 1 gX]) , 

where is the Euclidean action, r the imaginary time 
and the path integral is over all periodic trajectories in 
the interval [0, H0\. This path integral is Gaussian and 
can be evaluated exactly. It is convenient to represent 



Since a set of non-interacting oscillators can represent the most 
general Gaussian thermal reservoir, considering a non-diagonal 
mass matrix may appear superfluous. However we need the non- 



diagonal elements to generate the correlation function of two 
different reservoir momenta by partial differentiation of the par- 
tition function. 
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the integration paths via their Fourier series, which takes 
care of the condition on periodicity. 



x(r)= £ x r e— , 

I' — — OO 

oo 

X(r)= J2 X r e^ T , 



(C.4) 
(C.5) 



where x_ r = xj, X_ r = Xj (dagger stands for Her- 
mitian conjugation) since x(r) and X(r) are real and 
v v = Zirr/hfi are the bosonic Matsubara frequencies. 
Written in terms of the Fourier coefficients the Euclidean 
action becomes: 



S ( E) = M J2 (X r t(^ 2 M + C)X r ) + ^ ]T (xt, r 2 mx I+ [x 1 .-c- 1 gX r ] , c[x r -c- I gX r ] 



hp 

~2 



(C.6) 



Next we decompose x r = x r + y r where 



x r = (z^m + c) x gX r 



(C.7) 



is chosen such that does not have a term linear in 
y r . The action can be written as: 



c(E) = c(E) r 
reservoir lv 
oo 



hp 
hp 



? (E) 
^system 



;x], 



Y (yJ m + C ) Yr 



r= — oo 
oo 

, } (Xj( i / I 2 M + C + 2 i , 1 .7( i , r ))X r ) , 



r 



where the damping kernel is given by 



- / \ 1 T -i -1 

70) = 2 § m 2uJ 



(CIO) 



which is the Laplace transform of Eq. (II. 5). The parti- 
tion function of the closed system is given by: 

/oo 
n dX r exp(-^ ( y E s ) tom [X]/n) 

1--OG 

n. OO 

x J J] d yr exp(-^ c E s ) crvoir [y]/fi) . (C.ll) 



(C.8) The normalization factor Af is yet unspecified because it 
is not easy to determine the measure of the path integral. 
(C.9) J\f will be determined indirectly at the final stage of this 
calculation by considering the limiting case of no system- 
environment coupling. 

The integrals in Eq. (C.ll) are all Gaussian. Ignoring 
the normalization for now the integration gives: 



oo 

° * J}^ y/det \y? m + c] ^dct [v* M + C + 2 vfi&j] ' 

1 1 -pr 1 1 

" v/detR V^etlC] det m + c] det \y* M + C + 2 ^ 7 (i/ r ) 



(C.12) 
(C.13) 



In the second line we used the fact that the elements of 
the product corresponding to positive and negative val- 
ues of r are identical to restrict the product to positive 
r and pulled out the r — entry. To determine the nor- 
malization let us recall the partition function for a simple 



r 



harmonic oscillator: 



1 1 

Zm ° = 2sinh(/3ftu;/2) = Jhw JJ ^ + v 2 ■ ( C - 14 ) 
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This naturally generalizes to N harmonic oscillators by: 
1 1 



Z 



NHO = 



det[2sinh(/3/iu;/2)] /?ftdet[w] det[u> 



n 



(C.15) 



Using the definition (11.16) the partition function can also 
be written as: 

— J[|det(M^G(^)) . (C.17) 



In the limit of no coupling we demand that the partition 
function be a product of two partition functions of this 
form. This condition fixes the normalization and the final 
answer is: 



Zq = Zb x det 



1 



(C.16) 



OO / 

fldet 



r=l 



VO 2 + v$ + 2M-5y r 7(j/ r )M-5 



where Z-& — Tr [exp(— /3He)] is the partition function 
of reservoir oscillators without coupling to the system. 



Appendix D: Derivation of Eqs. (II. 33-11.38) 

In this appendix we derive some of the results pre- 
sented in Sec. II C. Angular bracket with the subscript 
C denotes expectation values in the closed system ther- 
mal state. Expectation values in the uncorrelated state 
are denoted by attaching the subscript E to the bracket. 
Note that the damping kernel depends on the environ- 
mental variables and the coupling constants alone. There 
is no dependence on system variables. Using Eq. (11.23) 
we calculate the single-time system position-position cor- 
relation as: 



<(xx T ) AB ) c 



1 d 

~PdC AB 



2 d 
PdC AB 



TrlogC+ ^T^^^^Trlog M^G^y 1 



I 9 



oo 



(m-.g ( ,,-.)-'m-.% 



1 2 - 

= ^G(i^o)ab + ^ 2J G (^)ab , 

(XX T ) c = iG(^ ) + |f:G(, r ), 



r=l 



(D.l) 
(D.2) 
(D.3) 
(D.4) 



where we used the fact that C and G(v r ) are symmet- momentum-momentum correlations can be calculated in 
ric matrices and G(u ) = G(0) = C _1 . The system a similar way using Eq. (11.25). 



<(PP T ) AB )c = 



' " Trlog^M- 1 ) + | g (M a i)AB E ( Trl °S t M_1 ] + Trl °S [G(^r)- 1 



M AB , 2 



M AB + Tr 



G(i/ r ) 



r=l 
^(M-^AB 



(p p t ) c = M + \ ( M - M M ) ■ 



(3 



r=l 



(D.5) 
(D.6) 
(D.7) 



We used Theorem 3 in the first line. In the second line For the system-environment position correlations note 
we used Theorem 4 for all terms and Theorem 5 for the that only the damping kernel depends on the interaction 
last term with A 1; A 2 ,A 3 M _1 , z/ 2 G(f r ), 1. 
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matrix: 
1 d log Z c 



1 OO 



9 



/3 5(g T ) A a /3^9(g T )Aa 

2 



TrlogCG^)- 1 ] , 



r=l BC 



Aa 



(D.8) 



The partial derivative of the damping kernel can be cal- 
culated explicitly. For this differentiation it is useful to 



rewrite ■y(i'i) as: 



2iyy(i^) = 2/ r 2 g T c _1 (m" 1 + c" 1 ) 1 c^g (D.9) 



For brevity of notation we define a(z/ r ) such that 7(V r ) 
§g T i/ r a(z/ r ) g. 



2 ^ ^ ^ = ^ ^ (g ) Cc (i/ r a(i/ r )) cf g£B , 



<9(g T ) A a <9gaA 



ef 



cf 



= 1 S ^ 5 ca (^ 2 a(z/ r )) cf g fB + (g T ) Cc (^ r 2 a(^ r )) cf 5fJ B A 



r a 0r) g) aB <*CA + (g T l' 2 a(zy r )) Ca £ B A 



(D.10) 

(D.ll) 
(D.12) 



Plugging this result in Eq. (D.8) we get: 
ld\ogZ c 2^ 

Using this result in Eq. (11.26) we get Eq. (11.35). 
To derive Eq. (11.37) we start from Eq. (11.32): 



2 d log Z c _ / nr jT\ 
~P dm 



r = (pp T ) 5 



(D.14) 



B ' am 1 

r— 1 



M G (l/ r ) 



I 

Using Theorem 4 we get 

a 



dim- 1 ) 



j— Trlog 



ab 



M~ 1 G _1 (i' r 



where 



Tr 



G(i/ r ) 



9G-^j = 9 7 (^ r ) 
9(m- 1 ) ab r 9(m- 1 ) ab 



aG- 1 ^) 

5(m- 1 )»b 
(D.15) 



(D.16) 



Next use Theorem 5 and plug the result back into 
Eq. (D.15): 



5(m ijab A<™ x > ■ 



Tr 



G(i/ r 



) 9G I V) 
9(m- 1 ) ab 



ABcd 



d(m 1 + tfc 1 
9(m -x )ab 

£ GK)ab ( -. 2 (g T a(, r ) c) fic 9(m ^ 1 +^ C ' 1)cd (c a, „ i gi 1A 



(rn-^^c- 1 )"^-^, 



ab 



(D.17) 
(D.18) 

(D.19) 



Observe that: 



It follows that 



9(m- 1 )ab 



Tr 



db ■ 



(D.20) 



G(z/ r ) 



aG- 1 ^) 



^(m-^ab 

= -E^ 2 G(^) AB (g T a(i/ r )c) Ba (ca(^.)g) 



bA ' 



AB 



= - (ca(i/ r )gi<G(i/ r )g a(z/ r )c 



ba. 



(D.21) 
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In the last step we used the fact that both G(y r ) and 
a(?/ r ) are symmetric matrices. Eq. (D.15) becomes: 







— Trlog M G (z/ r ) 



-ca(i/ r ) g f r 2 G(i/ r )g T a(^ r ) c . 



(D.22) 



We plug this into Eq. (D.14) and note that ca(i/ r ) = 
mf(i/ r ) to get Eq. (11.37). 

The derivation of Eq. (11.38) is almost identical to that 
of Eq. (11.37) but with more terms. We do not show the 
details of that derivation here. 



Appendix E: Proof of conclusions of Sec. II E 

Using the fact that all position-momentum correlations 
vanish we get: 

Xj) c = g T f (s) m (x Xj) c , (E.l) 
(£( S )Pj) c = g T f( S )( Po pT) c , (E.2) 

where the expectation values on the RHS are given by 
Eqs. (11.35,11.36). 



(|( S )XJ) C =g T m^ C °^m^g(X Xj) c (E.:>) 

,2 



m 



_i cos(u;s) v 2 
ijj 2 (u) 2 + v 2 



m 2gG(i/ r ). 



The hrst term on the right-hand side can be seen to decay 
by the fact that 



T _icos(o;s) 



m 2 



m 2 g = 2-y(s) 



(E.4) 



The second term can be seen to decay by noting the 
inequality 

T _i cos(u;s) v 2 _i T _icos(o;s) _i 

g m 2 2 , 3 , — jrm ^g^g^ 2 — m 2 g, 

u;^(u; 2 + f~) u> z 

< 27(5) , (E.5) 

in the sense of positive-definite matrix kernels, since both 
u: and (uj 2 + v 2 ) are positive matrices and cosine is 
a positive-definite kernel. The summation over r in 
Eq. (E.3) is finite as can be seen from Eq. (11.33). As 
a result (£(s)Xq) c is a function that decays over time 
like -f(s). When we take the convolution of this with 



another decaying function G(t — s) and let t — > 00 the 
overlap goes to zero. This way we argue that second line 
of Eq. (11.42) vanishes. A similar calculation establishes 
the same goes for the third line. 



= g T f (s) c (x Xj) M + g T f (s) g <X Xj) M , 



(E.6) 



1 OO 

r— 1 



s c m 2 



r=l 

1 X ^ d 

r— 1 



m 2 



o) 2 (a; 2 + v 2 ) 
sin(a;s) v 2 



m-3gG(z/ r )M. 



ljj(<jj 2 + ^ 2 ) 



(E.7) 

m-4 g G(i/ r )M, (E.8) 



1 cos ws vz 



m 2 



a; 2 (a; 2 + z/ 2 ) 



m 2 i 



G(i/ r )M. 
(E.9) 



The term inside square brackets decays as 7(5) as can be 
seen from Eq. (E.5) and the argument following it. The 
summation over r is finite as before. Hence (£(s)Pj) c 
decays over time like 'j(s). The convolution of this with 
another decaying function G(t— r) gives zero in the limit 
t — > 00. 

The second and third lines of Eq. (11.42) are zero for the 
uncorrelated initial state as well. This follows trivially 
from: <£(«) Xj) E = <£( a ) pT) E = 0. 

Finally we need to show that the fourth line of 
Eq. (11.42) is the same for both cases. This requires 
showing that the late-time limit of the noise kernel is 
the same. We know that the noise kernel is stationary 
for the uncorrelated initial state. Let us focus on the 
noise kernel of the closed system thermal state. 



(E.10) 



g T (f (a) m <x xj) c mf (a') + f (s) (poP?> c f (a 7 )) g • 



We use Eqs. (11.37,11.38) on the RHS. The derivation is 
straightforward but tedious. The theorems in App. B are 
utilized repeatedly. 

The uncorrelated noise kernel is obtained if only the 
first terms in Eqs. (11.37,11.38) are kept and the rest ig- 
nored. Hence we need to show that all the other terms 
vanish in the late-time limit. The strategy is the same as 
before: we show that these terms are bounded by a func- 
tion proportional to the damping kernel or its derivatives. 
We work out the details for two terms explicitly. 

First consider the term in the noise kernel Eq. (E.10) 
due to the second term in Eq. (11.38). 
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g T f (a) m c^g (X XT) c gTc^mf (a') g 

= g T m~ 2 cos(u;a) m _2 mm 2 o; 2 m 2 g (XqXq ) c g T m 2 o; 2 m 2 mm 2 cos(o>s') m 2 g . 



t _icos(u;sJ _i , T , T _i 
^m 2 2 g( X oX ) c g m 2 



s(us') 



m 2 1 



4 7 (a) (X XT) c7 (a') 



(E.11) 
(E.12) 
(E.13) 



Unlike previous cases we were able to express this term 
exactly in terms of the damping kernel. It is a decaying 
function in both a and a' variables. The convolution of 
7(a) with G(t — a) in Eq. (11.42) goes to zero if we let 



t —> oo. Similarly the overlap of -f(s') with G(i' — a') 
vanishes in the limit t' — > oo. 

Secondly consider the term in the noise kernel 
Eq. (E.10) due to the third term in Eq. (11.38). 



g f(a) mm 



m 2 gi/ 2 G(j/ r )g T m 2 



u 2 (uj 2 + v\ 



r u m 2 m 



f(»')g 



m 2 



cosfwa) i , T i cos (ws'l 

m 2 gz/ I 2 G(i/ r )g i m 2 „\ | 



M r=i 



u> 2 (u> 2 + V 2 ) 

T _i cos(u;a) _i 
+ ^) 



G(t/ r ) d 2 
da' 2 



m 2 g, 



T _i cos(u;a') i/ 2 _i 
S m 2 2 ^ 2 i ^ m 5 i 



(E.14) 
(E.15) 



As before we conclude that the terms in square brackets 
decay like the damping kernel. The summation over r 
is finite as can be seen from Eq. (11.33) and noting that 
v r > 1 for all positive r. 

Close inspection of all the other terms in Eq. (E.10) 
reveals that they have roughly the same form as those we 
worked out the details explicitly. All these terms vanish 
in the late-time limit. 

This proves the equivalence of the late-time limit of the 
uncorrelated initial state to that of the late-time limit of 
the closed system thermal state. Since the closed system 



thermal state is stationary our proof is complete. 
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